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We present numerically determined travelling- wave solutions for pressure-driven flow through a straight duct 
with a square cross-section. This family of solutions represents typical coherent structures (a staggered array 
of counter-rotating streamwise vortices and an associated low-speed streak) on each wall. Their streamwise 
average flow in the cross-sectional plane corresponds to an eight vortex pattern much alike the secondary flow 
found in the turbulent regime. 



I. INTRODUCTION 

The flow through a straight duct with rectangular 
cross-section shares important features with both circu- 
lar pipe flow and plane channel flow. Analogous to pipe 
flow, laminar flow in a square duct is believed to be sta- 
ble with respect to all infinitesimal perturbations 25 . On 
the other hand, for aspect ratios above a critical value of 
A crit R3 3.2, duct flow is linearly unstable. In fact, ducts 
with large aspect ratio are often used in laboratory exper- 
iments as a substitute for plane channels 26 . Despite these 
similarities, the fact that the cross-section of rectangular 
ducts is neither axisymmetric nor spanwise-homogeneous 
leads to an interesting phenomenon: turbulent flow in this 
geometry exhibits mean secondary motion in the plane 
perpendicular to its axis 27 . This mean cross-flow motion 
is directed towards the corners in the vicinity of the diag- 
onals, and away from the wall near the bisectors, giving 
rise to a pattern with eight vortices. Although of small 
amplitude (few percent of the primary flow), the sec- 
ondary motion significantly deforms the latter, thereby 
inducing strong variations of local mean wall-shear stress 
along the edges 28 . 

The results of recent direct numerical simulations of 
turbulent square duct flow suggest that the occurrence of 
secondary motion is a statistical footprint of the prefer- 
ential location of coherent structures at certain positions 
along the edges 29 . This conclusion is based upon data ob- 
tained through a vortex eduction study applied to snap- 
shots of turbulent flow fields at low Reynolds numbers. In 
contrast to these "approximate" coherent structures, "ex- 
act" coherent structures (travelling- waves, time-periodic 
solutions) in duct flow have only received little attention 
in the literature. To our knowledge, the only study pre- 
senting non-linear travelling- wave solutions for the square 
duct is ref. 30. However, in contrast to the claim of the 
authors, the solution family presented therein is not rep- 
resentative of the turbulent flow states observed in ref. 29. 
More specifically, the travelling-waves of ref. 30 exhibit a 
pattern of four mean secondary flow vortices which can 
appear in two orientations (rotated by 7r/2). In this case 
what should be compared is on the one hand the sta- 
tistically averaged turbulent flow field and on the other 



hand the average between the two (equivalent) travelling- 
wave solutions rotated by ir/2 with respect to each other. 
Performing this averaging on the solution of ref. 30, how- 
ever, results in a vorticity pattern with a sign which is 
opposite to the one of the mean secondary flow found in 
turbulence 29 . 

The objective of the present work is to determine exact 
coherent structures which are potentially relevant to tur- 
bulent flow in a square duct. For this purpose we have 
computed different families of travelling- waves. In this 
article we present one particular solution family which 
entails an eight-vortex secondary flow with the same 
symmetries and sense of rotation as that found in time- 
averaged turbulent fields. This enables us to directly 
compare the properties of the exact coherent structures 
to statistical quantities measured in turbulent flow. We 
are therefore in a position to give solid theoretical sup- 
port to our earlier conjecture 29,31 by showing that near- 
wall coherent structures are indeed capable of generating 
secondary flow in a square duct. 



II. NUMERICAL METHOD 

In the following the coordinates x = (a;, y, z) are cho- 
sen such that x is aligned with the duct axis and y, z 
are oriented along the duct bisectors. From the duct 
half-width h, the maximum laminar flow velocity u max 
and the kinematic viscosity v we can form the base flow 
Reynolds number Re = u max h/v. We will also use the 
bulk flow Reynolds number Re\, — u^h/v (based on the 
bulk velocity ttf,) as well as the friction-velocity based 
Reynolds number Re T = u T hjv (where the wall- friction 
velocity u T is defined in the usual way). The super- 
script + will be used to represent quantities normalized 
with v and u T . Introducing the decomposition of the 
velocity field u = U(y,z)e x + u'(x, t) and the pressure 
p = P(x) + p'(x,t), the flow equations for the perturba- 
tions read: 

d t u' + Ud x u! + (u' • V){Ue x + u') + Vp 1 = tA7 2 u', 

V • u' = 0, u'(x, y — ±h, z, t) = u'(x, y,z = ±h, t) = 0, 

(1) 



2 



(a) 

Nx c f £ 

2 0.0361062 8.5279e-3 
4 0.0364307 3.8270e-4 
6 0.0364184 4.5376e-5 
8 0.0364168 1.4451e-6 
10 0.0364167 5.6946e-8 
14 0.0364167 - 

(6) 

Ny = Nz Cf £ 

20 0.0320227 0.00526353 

30 0.0321351 0.00176935 

34 0.0321557 0.00113064 

40 0.0321759 0.00050299 

44 0.0321851 0.00021773 

48 0.0321921 - 



TABLE I. Convergence of the solution with the truncation 
level, (a) Varying the number of Fourier modes for fixed N y = 
N z = 30 {Re = 2345, ah = 1, lower branch), (fe) Varying the 
number of modified Chebyshev polynomials for fixed N x = 2 
{Re — 2600, ah = 1.5547, lower branch). The relative error is 
computed with respect to the result obtained at the highest 
truncation level c/ ire /, viz. £ = (c/ — c /. re /)/c/, r e/- 



where e x is the unit vector in the ^-direction, u' = 
(vf ,v' ,w'), the mass density has been taken as unity 
and dt = d/dt as well as d x — d/dx. The streamwise 
laminar base flow component U is solution of {d 2 jdy 2 + 
d 2 /dz 2 )U — Re ■ dP/dx with dP/dx = const. We em- 
ploy a primitive variable formulation, expanding the de- 
pendent variables ip — {u' , v' , w' ,p'} as follows (assuming 
the solution to be a travelling- wave) : 

N x 2N y -k^~> 2N z -k { ^ 

¥>(X.*)= J2 J2 Yl Vnml* 

n=-N x m=0 1=0 v ' 

x^\y)4> { *\z) exp(ina(x - ct)) , 

where a is the streamwise wavenumber, c the (real- 
valued) phase speed and i = The functions 
are modified Chebyshev polynomials which — in the case 
of tp being a velocity component — incorporate the wall 
boundary conditions 32 . No boundary conditions are pre- 
scribed for the pressure field, and hence two polyno- 
mial degrees less than for the velocity field are used, 
i.e. fc( p ') = 2 and fc("') = fc^') = = 0. The 
equations (1) admit solutions which exhibit particular 
combinations of odd/even cross-sectional parities as de- 
scribed below. By restricting our attention to these sym- 
metric solutions, the number of modes in the y and z 
directions can be halved without loss of resolution. 
After substituting (2) into (1), a Galerkin method is 
employed in the streamwise (Fourier) direction, and a 
collocation method is used in the cross-stream (Cheby- 
shev) directions. In our case, the collocation points for 
pressure are chosen as the Gauss points, i.e. pressure is 
staggered with respect to the usual Gauss-Lobatto grid 
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FIG. 1. Convergence of the solution curves for two different 
streamwise wavenumbers ah — {0.6, 1.5547} computed at dif- 
ferent truncation levels, (a) For ah = 1.5547, the dashed line 
corresponds to 2 x 30 x 30 and the filled circles to 4 x 32 x 32. 
For ah = 0.6, the truncation level varies along the arrow in 
the order: 2 x 30 x 30, 4 x 20 x 20, 6 x 26 x 26, 10 x 20 x 20, 
14 x 20 x 20, 20 x 20 x 20, with the filled squares indicating 
individual solutions at 20 x 26 x 26. (b) Close-up of the data 
shown in (a) for ah = 0.6 on the upper branch. The curve for 
the truncation level 20 x 20 x 20 is drawn up to the point where 
the relative difference with respect to the solution obtained 
at 14 x 20 x 20 is still less than one percent. 



u v w p 

I (o,e) (e,e) (o,o) (o,e) 

II (o,o) (e,o) (o,e) (o,o) 

III (e,e) (o,e) (e,o) (e,e) 

IV (e,o) (o,o) (e,e) (e,o) 

TABLE II. The four combinations of cross-sectional solu- 
tion parities admitted by the linearized operator, as shown 
in ref. 25. The notation is such that e.g. "(o,e)" stands for an 
odd parity in y and an even parity in z. 
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FIG. 2. Minimum wavelength (in wall units) of the travelling- 
waves given as a function of the bulk Reynolds number. The 
open symbols correspond to a truncation level of 2 x 30 x 30, 
the filled ones to 4 x 20 x 20. 



FIG. 4. Wall friction factor c/ = %u 2 T ju\ as a function of 
Reb in logarithmic scale (u r being the friction velocity). Line 
styles and symbols as defined in figure 3. The additional 
straight chain-dotted line indicates laminar friction. 



used for velocity, thereby avoiding a discretization of the 
continuity equation at the corner points. It should be 
noted that the above discretization of the space oper- 
ator guarantees a solution which is free from spurious 
pressure modes 33 . The solution of the resulting non- 
linear algebraic system is performed via Newton- Raphson 
iteration, the linear system being solved by means of 
LU decomposition, exploiting parallel computing tech- 
niques. Continuation along the different problem param- 
eters is achieved through an arc- length procedure 34 . It 
has been verified that our code reproduces the solution 
of ref. 30 with high accuracy (choosing ah = 0.85 leads 
to min(i?efc) = 596.1191 and c/u max = 0.4780 when us- 
ing N x = 4, N y = N z = 26). Furthermore, numeri- 
cal convergence of the solutions presented in the follow- 
ing has been thoroughly investigated. Table I shows the 
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FIG. 3. Three-dimensional perturbation velocity, normal- 
ized by the bulk velocity, as a function of the bulk Reynolds 
number. The lines connect solutions at constant streamwise 
wavenumber: ah — 0.6 (dashed), ah — 1 (chain-dotted), 
ah = 1.5547 (solid). The symbols (•) are for turbulent flow 29 . 



convergence with respect to the number of streamwise 
Fourier modes as well as the cross-stream polynomials 
for two particular parameter points. It was found that 
the truncation level required for resolving the present so- 
lution family varies considerably in parameter space. In 
particular, resolving the upper branch solution at small 
wavenumbers imposes the most stringent requirements 
(as previously observed in Couette flow 35 ). This point is 
illustrated in figure 1 , where solution curves at two differ- 
ent streamwise wavenumbers are depicted as computed 
with different truncation levels. The plotted quantity is 
the perturbation energy which excludes the streamwise- 
constant modes, i.e. E 3D — J2 n ^a 1 1 \u n \ 2 dydz where 
u„ is the nth Fourier coefficient of velocity. The figure 
shows that N x = 2 suffices for ah = 1.5547, while the 
smaller wavenumber value ah = 0.6 requires a consid- 
erably higher number of streamwise modes. In particu- 
lar, we have used a streamwise truncation level of up to 
N x = 20 on the upper branch for ah — 0.6. With a res- 
olution of N x — 20 and N y = N z = 20 we obtain results 
which are within one percent of the results obtained for 
N x = 14 and N y = N z = 20 up to Re = 6565 (corre- 
sponding to Ret, = 1247.2). Beyond that point, results 
obtained with these two truncation levels exhibit sub- 
stantially larger differences, and our finest solution will 
therefore only be discussed up to Reb — 1247.2 in the fol- 
lowing. These tough resolution requirements stem from 
the fact that the flow field on the upper branch at smaller 
wavenumbers develops very strong streamwise gradients 
as can be seen in figure 5 below. 

In the absence of a "natural" primary bifurcation point 
from which a reasonable initial guess can be constructed, 
we resort to the homotopy method proposed by Waleffe 36 
based on the so-called self-sustaining process of wall- 
bounded shear flows. Here streamwise vortices are ar- 
tificially added to the base flow and forced against vis- 
cous decay, leading to streaks, which are in turn linearly 
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unstable, feeding back into the original vortices. The 
non-linear solution is then continued back to the orig- 
inal problem, i.e. the unforced flow. As suggested by 
Waleffe 36 we select the streamwise rolls amongst the least 
decaying eigenmodes of the Stokes operator in the cross- 
sectional plane. In the following we will focus upon re- 
sults obtained using the fourth eigenmode which ex- 
hibits a parity of type III in the nomenclature of ref. 25 
(cf. table II). In fact, the non-linear equations admit so- 
lutions which are such that all even streamwise Fourier 
modes verify the parity type III and all odd Fourier 
modes simultaneously verify one of the parity types I, 
II, III or IV. This property allows us to reduce the sys- 
tem size while choosing any linear streak instability mode 
in conjunction with rolls of parity III. The following re- 
sults were obtained when starting at a primary bifurca- 
tion point of a linear perturbation with parity II. 



III. RESULTS 

The critical point (minimum value of the bulk 
Reynolds number) is found at min(i?eb) = 471 for 
ah = 0.927 (determined at a truncation level of N x = 4, 
N y = N z = 26); this corresponds to a value of the 
friction-velocity based Reynolds number of Re T = 51.6. 
It should be recalled that self-sustained turbulence in 
square duct flow exists above Ret « 1070 29 (Re T « 80); 
localized turbulent "puffs" have been found to subsist 
above Reb ~ 800 (unpublished results) . The region in pa- 
rameter space which is occupied by this travelling-wave 
family has an upper bound in wavenumber a which in- 
creases with the Reynolds number; when normalized in 
wall units the corresponding short wavelength limit A+ in 
is relatively constant with a value between 250 and 300 
over a considerable range of Reb (figure 2). This shortest 
streamwise wavelength at which the travelling-waves ap- 
pear is not very far from the length of the minimal flow 
unit for self-sustained turbulence which measures roughly 
190 wall units in the square duct 29 and 250-300 in plane 
channel flow 38 . It should be noted that a precise deter- 
mination of the upper limit of the streamwise wavelength 
as function of Reynolds number is beyond the scope of 
the present paper. 

Figure 3 characterizes the solution in terms of the 
three-dimensional perturbation energy. Here and in the 
following we will focus upon three specific values of the 
streamwise wavenumber, namely ah = {0.6, 1.0, 1.5547}. 

1 /2 

It can be seen that the quantity E 3D jub varies signif- 
icantly with the streamwise wavenumber (especially on 
the upper branch), the highest values being observed for 
the smallest of the three wavenumbers (ah — 0.6). The 
figure also includes the time-averaged counterpart of the 
three-dimensional perturbation energy for turbulent flow 
(obtained through DNS 29 ). At low Reynolds number 
the present solutions with wavenumber ah = {0.6, 1} are 
more energetic than the average turbulent flow, while the 
opposite is true for ah = 1 at higher Reynolds number 



Reb ^ 1370; unfortunately, due to the above mentioned 
high resolution requirements, it was not possible to fur- 
ther continue the curve for ah = 0.6 in order to judge 
whether a similar drop-off occurs at that wavenumber. 

Let us now turn to the friction factor. Figure 4 demon- 
strates that the lower branch of the present solution sur- 
face scales approximately as the laminar solution (i.e. 
varying as Re^ 1 ) at all wavenumbers. The upper branch 
exhibits significantly higher friction values at smaller 
wavenumbers. On the upper branch curves correspond- 
ing to ah = 0.6 and ah = 1 the average turbulent fric- 
tion values are exceeded for small Reynolds numbers; 
again, cross-over between the travelling-waves at ah = 1 
and the turbulent data takes place at Reb « 1370, be- 
yond which the upper branch solution scales approxi- 
mately as Re^ 1 . Conversely, for the larger wavenumber 
(ah = 1.5547) the upper branch has a similar slope to 
the lower branch at all Reynolds numbers and does not 
correspond to significantly increased friction with respect 
to the laminar base flow. 

In figure 5 solutions on the upper and lower branch at 
the same wavelength (ah = 1) and similar bulk Reynolds 
number (Reb ~ 1400) are depicted by means of isosur- 
faces of streamwise vorticity uj x and total streamwise 
velocity u. The former isosurfaces take the shape of 
streamwise elongated tubes which are slightly inclined 
with respect to the wall plane. Surfaces of the stream- 
wise velocity exhibit a single corrugation per wall, which 
corresponds to a low-speed streak. The streamwise vor- 
tices are in streamwise alternating (staggered) arrange- 
ment with respect to the nearest streak. When focusing 
upon structures in the vicinity of a single wall, it be- 
comes clear that the present travelling-waves are similar 
to exact coherent structures found in plane wall-bounded 
shear flows 36 . Furthermore, the staggered vortex/low- 
speed streak configuration described above is at the cen- 
ter of a multitude of observations in the near- wall region 
of turbulent flow (e.g. DNS of Jeong et al. 39 ) as well as 
in time-periodic solutions . 

The lower branch solution in figure 5 exhibits a very 
mild streamwise dependency (i.e. the streaks are nearly 
straight), as well as very smoothly shaped, large-scale 
streamwise vortices. This is consistent with the fact that 
the shape of the lower branch solution changes very little 
with Reynolds number. The upper branch, on the other 
hand, is characterized by very localized vortices and rel- 
atively sharp bending of the streaks. The high intensity 
of the vortices in the upper branch solution as well as 
their proximity to the walls is manifested through the 
existence of patches with strong mirror vorticity. 

Figure 5 also shows the streamwise-averaged flow cor- 
responding to the upper/lower branch solutions. A sec- 
ondary flow pattern with eight vortices, very much alike 
the one exhibited on average by turbulent duct flow, can 
be observed. The resulting deformation of the primary 
flow contours leads to a variation of wall-shear along the 
edges (not shown), which is characterized by two local 
maxima near the corners and a local minimum at the bi- 
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FIG. 5. Top row: the upper branch solution at Ret = 1404.1; bottom row: the lower branch solution at Reb = 1398.3; 
both at ah = 1. The graphs in the left and center column show surfaces of constant values of the total streamwise velocity 
u — 0.55 max(u) (sheet-like structure near the walls, green color) and of the streamwise vorticity at ±0.65 (±0.55) times its 
maximum absolute value in the case of the upper (lower) branch solution (tubular structures, blue/red color). This value 
corresponds to lj^ = ±0.0655 (wj = ±0.040). In the left column the surfaces are only shown on one side of a diagonal of the 
duct cross-section. The right column depicts the corresponding mean flow, the primary flow («) being visualized by contourlincs 
at {0.25,0.5,0.75} times the maximum value, and the secondary flow ({v), (w)) in form of arrows. Truncation level: N x — 20, 
N y = N z — 26. 
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sector - as found in turbulent flow 31 at Reynolds numbers 
below Rtb ~ 2000 . Figure 6 shows the statistically aver- 
aged mean flow from DNS 31 at Ret = 1400, facilitating a 
direct comparison with the present results of figure 5 at 
practically the same Reynolds number. A strong similar- 
ity of primary flow contours and secondary flow pattern 
between DNS data and the upper branch travelling-wave 
can indeed be observed. 

In order to quantify the amplitude of the streamwise 
averaged secondary motion the following energy norm 
can be defined: E vw — J J{\vo\ 2 + \wo\ 2 )dydz. The 
square-root of this quantity, normalized by the bulk ve- 
locity, is shown in figure 7. Analogously to the perturba- 
tion energy (cf. figure 3) , the present travelling- waves (at 
small wavenumbers) by far surpass the average secondary 
flow intensity of turbulence at lower Reynolds numbers. 
Again a steep drop occurs for ah=l around Re b « 1370, 
after which the travelling-waves and turbulence exhibit 
secondary flow of similar magnitude (to within 50% of 

1 /2 

each other). While the turbulent data for E v ' w /ub ap- 
pears to approach a constant value at higher Reynolds 
numbers, no scaling of the present upper-branch solu- 
tions can be deduced from our data. On the other 
hand, the curves for the lower branch practically col- 
lapse at all wavenumbers and vary approximately as 
Eyw /v,b ~ -Rer 1 , implying a scaling with the viscous ve- 

1 /2 

locity, i.e. E vw ~ v/h; this scaling has previously been 
observed in Couette flow 41 and pipe flow 42 . 

The phase velocity (shown in figure 8) corresponding 
to the present lower-branch solutions tends towards val- 
ues in the range of c/ub — 1.66 — 1.75 (depending on 
the wavenumber), which is comparable to the asymp- 
totic value 1.76 of the lower-branch solution with mirror 
symmetry in pipe flow 42 . On the upper branch and for 
small wavenumbers, the values for the phase speed of the 
present solutions are in the range of c/ttj = 1.2 — 1.4. 
In terms of the friction velocity, the lowest phase-speed 
c min ~ 12 is obtained near the turning point, i.e. the 
junction of lower and upper branches (for both wavenum- 
bers 0.6 and 1). This value is comparable to the propaga- 
tion speed of near-wall coherent structures in turbulent 
channel flow, which is equal to approximately 10 times 
the friction velocity 43 . 



IV. CONCLUSION 

The present non-linear equilibrium solutions to the 
Navier-Stokes equations are characterized by coher- 
ent structures (staggered streamwise vortices and wavy 
streaks) commonly found in wall-bounded shear flows. 
Their streamwise average also corresponds to the pat- 
tern of the time-averaged velocity field observed in (lab- 
oratory and numerical) experiments of turbulent square 
duct flow. Therefore, they demonstrate (for the first time 
theoretically) that near-wall coherent structures are in- 
deed capable of generating physically correct secondary 



motion in duct flow (as previously conjectured 29,31 ). Fur- 
thermore, we have shown that the present family of 
travelling-waves yields integral flow quantities which are 
generally comparable to statistical data of turbulent flow. 
This implies that these equilibrium solutions are poten- 
tially relevant to the turbulent state. In particular, it 
was shown that the upper branch solution at an interme- 
diate streamwise wavenumber (ah — 1) yields secondary 
flow intensities which are comparable to turbulence data 
for Reb > 1370; furthermore, the upper-branch solution 
changes from a highly energetic and high wall friction 
state at lower Reynolds numbers to a low-intensity state 
for higher Reb, crossing the statistical energy and wall 
friction level of turbulence around Reb ~ 1370. Inter- 
estingly, this drop in perturbation energy and wall fric- 
tion coincides with the Reynolds number range where 
turbulent duct flow exhibits a transition in terms of the 
instantaneous flow pattern: below Reb ~ 1250 the pre- 
dominant state of turbulent flow is characterized by a 
vortex pair and a single low-speed streak adjacent to 
one pair of opposite walls (the so-called 4- vortex state), 
whereas this feature gradually disappears with increasing 
Reynolds number, giving way to an 8- vortex state (vortex 
pair/streak present near all four walls simultaneously) 29 . 
Therefore, it can be expected that the present solu- 
tion family is pertinent to turbulence-induced mean sec- 
ondary motion in square ducts at Reb > 1370. Further- 
more, DNS shows that the average number of streaks 
which are simultaneously present along one edge of the 
cross-section increases with the Reynolds number 31 , con- 
sistent with a scaling of the spanwise streak-spacing 
in wall units. This property implies that additional 
travelling-waves with a distinct vortex pattern (16, 24, 
32,... vortices in total, i.e. 5, 7, 9,... streaks per edge) are 
expected to come into play progressively as the Reynolds 
number is increased. As a consequence, we expect our 
present family of exact coherent structures to be consis- 
tent with actual buffer-layer structures in turbulent flow 
only around Reb ~ 1370. 



It should be mentioned that other families of travelling- 
waves which are also consistent with the coherent buffer- 
layer structures and the average secondary flow pattern 
encountered in square duct turbulence might very well 
exist. However, the central point of the present commu- 
nication is the fact that such a non-linear solution does 
indeed exist at all. 



In order to determine which role the present equilib- 
rium solutions play in the dynamics of the initial value 
problem, the stability characteristics of the travelling- 
waves should first be investigated. Then actual turbu- 
lent phase-space trajectories need to be mapped out, and 
their relation to the travelling waves (fixed points) needs 
to be established (as e.g. performed for Couette flow in 
ref. 44). This analysis will be the object of a future study. 
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FIG. 7. The secondary flow intensity of the present solution 
family, normalized by the bulk velocity and shown as a func- 
tion of Reb- Line styles and symbols as defined in figure 3. 
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FIG. 8. The phase velocity normalized by the bulk velocity. 
Line styles as defined in figure 3. 



